Kadanoff-Baym description of Hubbard clusters out of equilibrium: performance of 
many-body schemes, correlation-induced damping and multiple steady states. 
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We present in detail a method we recently introduced (Phys. Rev. Lett. 103, 176404 (2009)) 
to describe finite systems in and out of equilibrium, where the evolution in time is performed via 
the Kadanoff-Baym equations within Many-Body Perturbation theory. The main non-equilibrium 
property we analyze is the time-dependent charge density. An other quantity we study is the 
exchange-correlation potential of time-dependent Density Functional Theory, obtained via reverse 
engineering from the time-dependent density. Our systems consist of small, strongly correlated 
clusters, described by a Hubbard Hamiltonian within the Hartree-Fock, second Born, GW and T- 
matrix approximations. We compare the results from the Kadanoff-Baym dynamics to those from 
exact numerical solutions. The outcome of our comparisons is that, among the many-body schemes 
considered, the T-matrix approximation is overall superior at all electron densities. Such compar- 
isons permit a general assessment of the whole idea of applying Many-Body Perturbation Theory, 
in the Kadanoff-Baym sense, to finite systems. A striking outcome of our analysis is that when 
the system evolves under a strong external field, the Kadanoff-Baym equations develop a steady- 
state solution as a consequence of a correlation-induced damping. This damping is present both 
in isolated (finite) systems, where it is purely artificial, as well as in clusters contacted to (infi- 
nite) macroscopic leads. To illustrate this point we present selected results for a system coupled 
to contacts within the T-matrix and second Born approximation. This damping behavior is in- 
trinsically linked to the Kadanoff-Baym time-evolution and is not a simple consequence of possible 
limitations/approximations in the calculation of the initial state. The extensive numerical charac- 
terization we performed indicates that this behavior is present whenever approximate self energies, 
based upon infinite partial summations, are used. A second important result is that, for isolated 
clusters, the steady state reached is not unique but depends on how one switches on the external 
field. Under some circumstances this is also true for clusters connected to macroscopic leads. We 
provide some statements of more general and conceptual character on how the damping mechanism 
depends on the system size and the number of particles, and conclude with an outlook and glimpses 
of future work. 

PACS numbers: 71.10.-w 71.10.Fd 73.63.-b 
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I. INTRODUCTION 

The Kadanoff-Baym equations (KBE)i— are one of 
the fundamental theoretical schemes of a microscopic de- 
scription of quantum systems out of equilibriumi^i^ Due 
the growing interest in time-dependent phenomena, in re- 
cent years the KBE have been the object of considerable 
attention in several branches of physics. Notable applica- 
tions of the KBE are in the areas of molecular quantum 
transport, high energy coupled plasmas, nuclear matter, 
astrophysics, to mention a fewi^"— Another favorable el- 
ement to the widespread use of KBE is the constantly ex- 
panding capability of today's computers, which has made 
the full numerical solution of the KBE possible. A main 
strength of the KBE is that one can, in a constructive 
way, build approximations of increasing complexity for 
the one-particle Green's function, G, the key quantity in 
the KBE. These approximations are obtained via Many- 
Body Perturbation Theory (MBPT) and are known as 
conserving approximations, since they guarantee the con- 
servation of important quantities such as total energy, 
number of particles, linear and angular momentum. 

However, the fulfillment of such conserving conditions 
is no guarantee of the quality of the actual results ob- 



tained within a specific MBA. Hence, it would be useful 
to have a way to assess the performance of a given con- 
serving MBA. One of the aims of this paper is to evaluate 
the range of validity of a group of well known Many-Body 
Approximations (MBA:s) by comparing the one-particle 
densities with the exact results for finite strongly corre- 
lated clusters out of equilibrium. The main attractive 
feature of such comparisons to exact results is the possi- 
bility of scrutinizing the performance of the Many-Body 
Approximations in the non-equilibrium regime. This 
knowledge is most valuable if one wishes to use approxi- 
mate many-body schemes for systems (typically, infinite 
ones, e.g. systems coupled to macroscopic leads as those 
we consider towards the end of the paper) where exact 
solutions are not available. 

Describing small clusters with strong electronic corre- 
lations, and subject to time-dependent fields, presents in- 
terest not only from a conceptual point of view; there are 
very many instances where the technological relevance of 
cluster physics is manifestly evidenti ^^d^ The main focus 
of this work is to investigate the dynamics of finite clus- 
ters. In this regard, we provide here a detailed account of 
a study performed very recently,— producing additional 
results for clusters and presenting in great detail the 
methodology we developed. We will, however, present 
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some results for a system contacted to macroscopic con- 
tacts to generalize some of our main findings. Similar 
clusters to those discussed here, coupled to leads have 
already been considered either in the stationary limil^ii 
or in the real-time domainJ^ii^ 

As specific finite model systems, we consider open- 
ended linear chains^ with short ranged Hubbard inter- 
actions. We study their dynamics through exact diag- 
onalization methods and by propagating the KBE for 
different MBA:s. The approximations we consider are 
the Hartree-Fock (HFA), the second Born (BA), the GW 
(GWA)^ and the T-matrix (TMA),^i^ Ah these ap- 
proximations are conserving,— which clearly is of great 
importance when propagating the KBE, and all of them, 
apart from HFA, have self energies which are non-local 
in space and time. For the GWA, we will consider both a 
spin-independent and a spin-dependent versionJi^^ The 
latter has the advantage to alleviate the effect of self- 
screeningi2ii2^ 

A general outcome of our study is that the TMA is seen 
to perform better than the other MBS:s at all fillings and 
interaction strengths. 

With the exact and KBE cluster dynamics at our 
disposal, we also investigate numerically a well es- 
tablished relation between MBPT and another frame- 
work to treat non-equilibrium phenomena, i.e. time- 
dependent Density Functional Theory (TDDFT)^i^ 
Using a spin-independent TDDFT description for the 
spin-compensated Hubbard model;^ we will obtain the 
exchange-correlation (xc) potentials corresponding to the 
different MBA:s via reverse engineering from the time- 
dependent densities. 

Our results show that the time-dependent KBE present 
two interesting features. The first is that for large exter- 
nal fields, the KBE time evolution in clusters exhibits 
a damped behavior induced by many-body correlations 
(hereafter, we refer to this as correlation-induced damp- 
ing). The second feature is that the steady state one 
reaches is not unique but depends on how the pertur- 
bation is switched on. We also investigate clusters con- 
nected to macroscopic leads, which we in this paper we 
treat within the TMA and BA. In this case the presence 
of the correlation-induced damping is in general not ar- 
tificial but the steady states may be non-unique. 

These statements are also valid in presence of macro- 
scopic contacts, as shown below for a small cluster cou- 
pled to leads and with a time evolution within the TMA. 

In finite clusters the correlation-induced damping and 
the existence of multiple steady states is artificial and we 
show that it is due to limitations of self-consistent Many- 
Body Perturbation theory when applied to finite systems. 
We will also provide selected numerical examples, which 
suggest that correlation-induced damping can be present 
even in the presence of leads. 

The question whether or not the correlation-induced 
damping and the existence of multiple steady states in 
contacted systems is a mere consequence of MBPT is at 
present not so straightforward to address with full gen- 



erality. 

The paper is organized as follows: we start with a de- 
scription of our model system(s) in section [III then, in 
section lHll we discuss the general properties of the single- 
particle Green's function. Section IIVI is an overview of 
the Many-Body Approximations used in this work. Sec- 
tion |V] is devoted to the procedure to obtain the ground 
state within the KBE. How to solve the KBE for the time 
evolution is reported in section IVIl In section I VIII we 
detail how we extract the TDDFT exchange-correlation 
potentials corresponding to a chosen MBA. The ground- 
state and time-dependent results are presented in sec- 
tion IVIIII and IIXI respectively. Section |X] deals with 
the correlation-induced damping which occurs during the 
KBE time evolution and the existence of multiple steady 
states. Finally, in section IXll we present our conclusions 
and direction for possible future work. 



II. MODEL SYSTEMS 

We will consider one-dimentional clusters with AI sites 
and with one orbital at each site. Thus, each site (or, 
equivalently, each orbital) can accommodate a maximum 
of two electrons, with opposite spin. The clusters may 
either be isolated, in which case the Hamiltonian in stan- 
dard notation is (we set the on-site energies equal to 0) 

Hc^ -V 

{RR').(7 R R.o- 

(1) 

or they may be attached to non-interacting leads of infi- 
nite size, as discussed below. In Eq. ([T]), fiRcr = a\^^aRa, 
a =t,i, and {RR') denotes pairs of nearest neighbor 
sites. The hopping parameter V = 1 and wr, (i) is a 
local external field which can be of any shape in time t 
and space. U and wr (t) are given in units of V. We 
will consider clusters with M = 2, 4, 6 sites and, without 
leads, Ne = 2, 6 electrons (in the presence of leads, the 
average number of electrons in the clusters is in general 
non- integer). Our approach is valid for systems which are 
compensated as well as uncompensated in spin. However, 
in what follows wc will only consider clusters (with/out 
leads) with an equal average number of spin-up and - 
down electrons in the ground state; this will hold at 
all times during the dynamics, since H has no spin-flip 
terms. Henceforth, n = = n^, where, n„ = N^/M 
and Ncr = (^r^c^R'^)- presence of leads, L, the 

Hamiltonian is 

H = Hc + Hl + Hlc (2) 

where describes non-interacting one-dimentional semi 
infinite chains, 

Hl = - Vl ^ al^^QR'a, (3) 

{RR'), a- 
R,R'eL 
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M) 



FIG. 1: Keldysh contours. 



and Hlc describes hopping between the central region, 
C, and the leads 



Hlc = -Vlc X! "flfT^fl'cr + h.c. 



(4) 



(RR'), a 
ReL,R' ec 



For isolated clusters, we use a short iterative Lanczos 
propagation to obtain the exact time evolution. A de- 
scription of our approach for approximate solutions in 
isolated and contacted clusters is the object of the next 
five sections. 



III. THE ONE-PARTICLE GREEN'S 
FUNCTIONS 

The one-particle Green's function is a reduced quan- 
tity, containing much less information than the underly- 
ing wave function. In general, it describes a system con- 
nected to a bath with which energy and particle can be 
exchanged. The knowledge of the Green's function gives 
access to the expectation values of all single-particle op- 
erators, excitation energies and the total energy of the 
system. 

The general definition of the one particle Green's func- 
tion is 

g {ri(TiZi,Y2(J2Z2) = (riCTil^ {zi,Z2) |r2Cr2> (5) 



i (U (-i/3,0)r^ iJjh (ricriZi)-0}^ (r2cr2^2) 



(6) 



where (...) denotes expectation values of the equilib- 
rium ensemble, ^pH and tplf are the field operators in 
the Heisenberg picture, 7 is the Keldysh contour, see 
Fig. mi,ii), is the path ordering operator, U is an 
evolution operator and f3 = l/ksT is the inverse tem- 
perature of the bath. The r and the s denote the la- 
bels corresponding to space (site) and spin coordinates of 
the single-particle basis. In this basis the Green's func- 
tion becomes a matrix. The variable z belongs to the 
Keldysh contour and is in general complex. For nota- 
tional convenience we denote real (imaginary) times by t 
{it) . From the definition of the Green's function one can 
derive the so-called Kubo-Martin-Schwinger condition;^ 
Q{zi,Z2) = -G {zi,Z2±il3). 

The expectation values of all the single-particle oper- 
ators are obtained according to 



(A it)) 



itt[a it)g{t,t+)] , 



(7) 



where the trace, Tr = J2rir2<7a' ^rir2^cTcr', is over space 
and spin indices. As discussed later and in the Ap- 
pendix A, the total energy can be found by evaluating 
the Galitskii-Migdal or Luttinger-Ward functionals. 

Since in this paper we study only paramagnetic sys- 
tems in spin-independent fields, one-particle quantities 
such as the G and the E become spin-diagonal: 



g (riCTiZi, r2(T2Z2) = G (ri^i, r222) <5aic 



(8) 



In the following, we will use the shorthand notation, 
1 = (ri,zi), etc., for space time coordinates. 

The Green's function obeys an integral equation (the 
so-called Dyson equation) 



G (12) = Go (12) + / Go (13) S (34) G (42) d34, (9) 
with the non-interacting Green's function Go defined by 



(z9,, -/i(l))Go(12) = <5(12) 



(10) 



In Eq. ((TU)) . h is the non-interacting Hamiltoniani^ It 
is convenient to decompose h as h{t) ^ t + w (t) ~ /i, 
with i ) i the one-particle kinetic energy (— Vf/2 in co- 
ordinate space) , ii) w (t) a local external field which may 
depend on time, and iii) fx the chemical potential. The 
latter is taken to be in between the last occupied and the 
first non occupied levels. The inclusion of ^ in /i implies 
that the Fermi energy is placed at zero energy; electron 
excitations thus have positive energies, while hole states 
have negative. The kernel of the Dyson equation, E, 
is called the self energy and is in general non-local in 
space and time. In the exact theory as well as in con- 
serving approximations^ the self energy is a functional 
of the Green's function, E [G] , and the Dyson equation 
must thus be solved sclf-consistently. In equilibrium all 
quantities depend only on z = Z2 — zi and the equations 
are then most easily handled, in terms of Fourier trans- 
formed quantities, in the frequency domain. There, the 
corresponding Dyson equation becomes a simple matrix 
equation, involving matrix multiplications; in the single- 
particle basis. 



G(e) = Go (e) + Go (e)E(e)G(e). 



(11) 



In this paper, we will work only with systems for which 
the initial state is the actual ground state (i.e. /3 —> 00). 
That is, thermal mixtures (/3 < 00) of initial states will 
not be considered. In equilibrium, all two-point propa- 
gators can be expressed in terms of a spectral function. 
Specializing to the Green's fimction, the spectral decom- 
position has the form 



G(e) 



A{e') 



e + 11] sgn (e') 



(12) 



where the spectral function, A (e) is related to the anti- 
Hermitian part of the corresponding propagator, which 
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for the G is: A (e) = -tt"! [G (e) - G'^ (e)] sgn (e). The 
Fermionic spectral function are positive definite and the 
one for G is normahzed: 



(13) 



where 1 represents the identity matrix in the single- 
particle basis. 



IV. MANY-BODY APPROXIMATIONS 

In general one can not construct the exact self energy 
and thus needs to rely on approximate schemes. In Many- 
Body Perturbation theory (MBPT), one can systemati- 
cally construct self energies of increasing complexity. The 
main idea is making a diagrammatic expansion of the self 
energy, and selecting different classes of diagrams which 
are then summed up to infinite order. There is a very 
important group of approximations which conserve quan- 
tities such as the total energy, the number of particles, 
linear and angular momentum, when the system is sub- 
ject to external fields. These conserving properties are 
related to the fact that the self energy is a functional 
derivative of a generating functional 



S(12) 



(5$ 



5G [21) 



(14) 



The use of conserving approximations is in general very 
important and, in fact, practically mandatory when 
studying non-equilibrium phenomena. In this paper we 
will study the conserving Hartree-Fock, second Born, 
GW and T-matrix approximations (HFA, BA, GWA 
and TMA respectively), see Fig. [2] The bare inter- 
action is taken to be local in time, Z-/ (ri, r2) (5 (ti, t2)- 
The formalism we will present is general but we remind 
the reader that in this paper we will only consider lo- 
cal interactions, Z^(ri,r2) = U5{rx,Y2,) When special- 
ized to our Hubbard clusters with one orbital/site (de- 
noted by R), the on-site interaction can be treated either 
as spin-dependent, U ^nnnfui or as spin-independent 

Y^Rao' °^Ra"^Rcr''^R'y"^R'y Thcsc two ways are evi- 
dently equivalent in any order by order expansion such as 
the HFA or the BA. In approximations based upon par- 
tial summations, however, this equivalence may be lost. 
To illustrate this point we consider the GWA both spin- 
independently (GWA) and spin-dependently (SGWA). 
The TMA is treated only spin-dependently. 

Hartree-Fock approximation. The simplest many-body 
treatment is given by the Hartree approximation (HA), 
where one takes only the first order direct term into ac- 
count (i.e. the exchange is excluded). Due to its rather 
simple nature it will not be consider further. 

The Hartree-Fock approximation (HFA) includes also 
the first order exchange diagram, the Fock term. The 
inclusion of this diagram, among other things, cures the 
self interaction of the HA. The resulting self energy is 



HA) 



HFA) 



BA) 



GWA) 



SGWA) 



TMA) 
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FIG. 2: Diagrammatic expansions of Many-Body Approxima- 
tions. 



local in time and thus constant in frequency space. The 
remaining diagrams are responsible for the many-body 
correlations and give rise to self energies which are non- 
local in time. 

It is often convenient to separate the time-local HFA 
and correlation contributions and write 



S (zi,Z2) = SffF 6 (zi, Z2) + (zi, 2:2) 



(15) 



Second Born approximation. The simplest scheme 
which involves correlations is the second Born approx- 
imation (BA) which corresponds to keeping all the dia- 
grams up to second order. 

GW approximation. The GW approximation is the 
leading term in the expansion of the self energy in terms 
of the dynamically screened interaction W . The expres- 
sion for the self energy in time space is given by 



Hgw (12) = Sff + (12) W (12) . 



(16) 



It should be noted that this expression does not involve 
any matrix multiplication. The screening of the bare in- 
teraction, U , results from all possible electron-hole ex- 
citations which are described by a series of bubble dia- 
grams, involving an irreducible polarization propagator, 
P. This series can be summed, yielding in frequency 
space, for a spin-independent interaction. 



W =U+UPW 
and, for a spin-dependent interaction, 
W = UPU + [UPf W. 



(17) 



(18) 



We remind the reader that these Dyson like equations 
involve matrix multiplications. In both cases^ the po- 
larization propagator, in time space, is 



P(12) = -iG(12)G(21). 



(19) 
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T-matrix approximation. The T-matrix approxima- 
tion comes from building the T-matrix, T, by summing 
all the ladder diagrams, representing electron-electron or 
hole-hole scatteringi^ The expression for the self energy 
is given by 

Stm {12)^T.HF+iJ Z^(13)G(43)r(34)Z^(42)d34. 

(20) 

In the case of an on-site, sitc-indcpendcnt interaction this 
simplifies to 

Etm (12) = ^HF + lU^G (21) T (12) . (21) 

The sum of the ladder terms in the T-matrix results in 

T ^(j)- <f)UT, (22) 

where the so-called irreducible vertex (j) is defined as 

</)(12) = -iG(12)G(12). (23) 

V. GROUND STATE 

The ground state is obtained by solving the Dyson 
equation, Eq. (|lip self-consistcntly. For clusters con- 
tacted to non-interacting leads, the problem can be ex- 
pressed entirely in terms of propagators which refer only 
to the central region and an embedding self energy j^^i^i 



Se™6(e) = ^|^LCp.9L (e) 



(24) 



where gL (e) is the non-interacting Green's function of 
the uncontacted lead L. The full Green's function in 
the central region will now obey a Dyson equation with 
both a many-body and an embedding self energy. The 
presence of Semb gives rise to continuous spectra, and 
standard techniques can be used to find self-consistent 
solutions. 

For the isolated clusters, we used a meromorphic rep- 
resentation to be described below. 



A. Meromorphic representation of finite systems 

For convenience all one-body quantities are rep- 
resented as matrices in a single-particle basis, e.g. 
Grri (e) = {R\G{e)\R'). In a finite system with a fi- 
nite phase space, all the spectral functions arc discrete: 



Arr. (e) = ^A]^^,5(e-a,), 



(25) 



where A-'j^j^, is a residue matrix and Oj is a pole position. 

^From Eq. we sec that the propagators themselves 
become meromorphic. 



Grr> (e) = E — 



A^RR' 



e — Uj + irj sgn (e) 



(26) 



One main advantage in using a meromorphic repre- 
sentation is that convolutions and cross-correlations are 
made analyticallyi^ Given the two functions 



Arr^ (6) ^ ^ , Brr, {e)^Y.— 



then their cross-correlation 



Crr.. (e) = j Arr, (e') Brr, (e + e') 



becomes 



Crr, {e) = Y,A'kR,B^ 



RR' 



jk 



e + Gk - bj 



(27) 



(28) 



(29) 



A second important attractive feature of a meromorphic 
representation is that one can compute at once the equi- 
librium many-body quantities with any time argument, 
both real and imaginary. For our Hubbard clusters, each 
of the quantities G, S, P, W^, $, T will be expressed in such 
representation during in the actual calculations. 



B. Solution to the Dyson equation 

The solution to the Dyson equation for the G can for- 
mally be written^ 

G (e) = [G„-i (6) - S (6)] -' = [e~h-j: {e)]-' . (30) 

The solution of this matrix equation is obtained in two 
steps^: 1) search of the pole position and 2) calculation 
of the residue matrix. 

1) The pole positions of G correspond to the zeros of 



det[e-h-Y. (e)] , 



(31) 



which we find by ordinary root finding algorithms. 
2) Once the pole positions are found we calculate the 
residue matrices by integration in the complex plain. We 
have in general 



^ / (e) = 27ri ^ Res (/, aj) , 



where 



/(^)-E 



A. 



(32) 



(33) 



If we now perform a closed integration around the pole 
Oj we obtain directly the residue matrix 

1 



27ri Ja, 

This integration is, in practice, performed numerically. 
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C. Self-consistency 

To reach self-consistency, we start by constructing the 
self energy with some initial G, normally taken to be Go, 
and then solve the corresponding Dyson equation. The 
resulting G is then used to build the now self energy and 
the procedure is carried on until convergence. In order 
to keep the number of poles under control (such number 
increases rapidly from iteration to iteration), we make 
use of a "decimation" procedure, where poles are merged 
if there arc small or close enough. When two poles are 
merged, the new pole position is the old center of mass 
position ("mass" being the trace of the residue matrix) 
and the new residue matrix is the sum of the two old 
residue matrices. To improve the convergence we update 
G by making a linear combination of the new (solution of 
the Dyson equation) and the old G:s. There are different 
functionals which yield the total energy of a system, all 
of which are equivalent at the point of stationary solu- 
tion to the Dyson equation. However, away from self- 
consistency, they do in general not coincide. One inde- 
pendent way of evaluating the degree of self-consistency 
is thus by comparing the values of different energy func- 
tionals, see appendix A. 



VI. TIME DEPENDENCE 

When an external field is applied to a system, it is in 
general driven out of equilibrium. When the system is out 
of equilibrium and away from the steady-state regime, all 
the quantities will intrinsically depend on the two time 
arguments (ti,i2) separately and the Keldysh formalism 
becomes essential. To explicitly show which convention 
we used in this paper for the path-ordered two-point 
Keldysh functions, /C, we give some definitions. 
A general Keldysh function can be written 

/C(12) = /C*(12)<5(zi,z2) 

+e(12)/C>(12) + e(21)/C<(12), (35) 

where > (<) refers to the hole (electron) part, and /C* 
the time-local part. 

The time coordinates are on the Keldysh contour and 
may be real or complex. For real times we may also 
introduce retarded and advanced propagators, 

/C^'(12) = /C^(12)5(zi,Z2) 

+e(ii,t2) [/C> (12) -/C< (12)] , (36) 

/C^(12) = /C^(12)(5(zi,z2) 

-e(i2,ti)[/C>(12)-/C<(12)]. (37) 

Note that for the /C:s considered in the paper, only E and 
W have parts which arc local in time. 



When both time arguments arc imaginary, the Keldysh 
function reduces to the corresponding equilibrium Mat- 
subara function: 

K.'^^ {T-T') = -lJCi-iT,-iT'). (38) 

In Eq. 9 (12) should be understood as 9 (zi, Z2), a 

generalized Heaviside function for zi , Z2 on the ordered 
Keldysh contour. In is worth noting that when both time 
arguments lie on the imaginary (Matsubara) axis, the 
quantities represent the initial, equilibrium, state which 
depend only on the time differences. As an example of 
how these terms are found in equilibrium in the mero- 
morphic representation we display the hole contribution 
to the Green's function: 

G< (ii,i2) XI ^J'^'"' (^3^) 

when both time arguments are real, 

G< {t, -ir) = i Aje'^^^e"'^ (40) 

when both one time argument is real and one imaginary, 

G< (-iTi, -iT2) = i X Aje-"^^^^-^^'^ (41) 

when both time arguments are imaginary. 

A. General symmetries 

We recall some prominent symmetry relations which 
will be used during the time propagation. ^^From the def- 
inition of the Green's function, Eq. ((5]) , one can derive a 
very important symmetry;^ which enters many relevant 
and useful relations: 

G^ (12) = -G^ (21)^ . (42) 

Additionally, from the definition of the retarded and ad- 
vanced Green's functions we obtain 

G^'^ (12) = G^'^ (21)^ . (43) 

From the expansion of the T and W in terms of and P 
it follows that the T and W will have the same symmetry 
properties as cf) and P. The symmetries of 4> and P can 
be deduced from their definitions 

P^ (12) = -P^ (21)^ =^ (12) = (21)^ (44) 

and 

0=§ (12) = -0^ (21)^ =^ T«; (12) = -T^ (21)^ . (45) 

In a similar way we find a symmetry relation, valid for 
all approximations, for the self energy, 
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(12) = (21)^ . (46) 
An additional symmetry fulfilled by W is 

W^{12) = W'^{21), (47) 

which implies 

W> (z, z) = W< (z, z) , Re Wrw (z, z) = V i?, i?'. 

(48) 

No equivalent relation exists for T . 

B. Solving for the non-equilibrium Green's 
function 

To obtain the non-equilibrium Green's function we 
need to solve the corresponding equations of motion, 
called the Kadanoff-Baym equations (KBE) , 

{idt, - h{l))G {12) =5 {12) + j S (13) G (32) d3, (49) 

{-idt^ - h {2)) G {12) ^5 {12) + ( G (13) S (32) d3.(50) 

The kernel of these equations, the S, will in general be 
a contraction of the Green's function with an other quan- 
tity which involves an infinite order summation such as 
the T or W . These quantities are defined by correspond- 
ing integral equations. Specializing to the case of T, we 
have 

T(12) = $(12)- /$ (13)^^(34) T (42) d34, (51) 
T (12) = $ (12) - /r (13) Z./ (34) $ (42) d34. (52) 

J 

We thus have two sets of coupled integral equations which 
need to be solved simultaneously at all times. 

C. Solution of coupled integral equations 

From the symmetry relations Eqs. (|^ we see 

that G and W, T are only needed on the upper/lower time 
matrix. We choose the lesser components on the upper 
triangle ti > t2 and the greater ones on the lower trian- 
gle ti <t2^ Propagation is thus made by expanding the 
Keldysh functions on the time square from T to T + A, 
see Fig. [3l In the integral equations, Eqs. (|49l [50l [51] 
I52p . the function to be determined appears in both the 
right and left-hand sides. To solve these equations we 
use a self-consistent predictor-corrector method.— The 
method can be described schematically by an external 
loop, where an approximate G at T -I- A is generated by 
propagating Eqs. (^51 [501) . ^'^'^ internal loop, where 



r + A T 




FIG. 3: (Color online) Time square. 

Eqs. (|5T1 [52)) are solved self consistently for a fixed ker- 
nel (j) [G] . The external loop is performed by a predictor- 
corrector method described below while the internal one 
is solved by the same iteration method as described for 
the ground state. The external loop is initiated by an ex- 
trapolated value of the collision integrals S^G, whilst 
the internal loop is started by an extrapolated value of 
the VF, T . The new time step is generated when the 
external loop achieves self-consistency. 



D. Kadanoff-Baym equations 

When propagating the Green's functions it is con- 
venient to separate the terms which are local in time 
and single-particle like {h and '^hf) from the remaining 
correlation-induced ones^ and introduce 

\}^h + i:HF- (53) 

In this way, f) subtitutes h and the full self energy is 
replaced by its correlation part Ec in the KBE (c/. Eqs. 
(UnilSni))- The reason for this partitioning is twofold: on 
the one hand the contribution from the single particle 
evolution is very important (and could thus lead to large 
numerical errors in the correlation contribution) and on 
the other it can be solved essentially in an exact way as 
will be detailed in Sec [VTfI 

There are several equivalent contours on which one can 
define the KBE. We use the contour ii) in Fig. |TJ which is 
numerically more stable and has an analytical limit when 
the temperature goes to zero. Once we specialize to this 
contour, the KBE become 

idt,G^{hM) ^l]{h)G^{ti,h)+lf{tuh), (54) 

-idt.G^ {h,h) = G^ {h,h) t) (<2) + if {tiM) ,(55) 
idtG< {t, -ir) = t) {t) G< (t, -ir) + I< (t, -ir) , (56) 
-idtG> {-iT, t) = G> {-iT, t) f) {t) + I> {-ir, t) (57) 

where the collision integrals with both time arguments 
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real are 



macroscopic leads, the problem can again be expressed in 
propagators which refer only to the central region, and an 
embedding self energy which now depends on timer^^ 



P/2 

'1 ' 



(ti , -ir) G^i-iT, ta) + (ii , ir) G^ [if, is) 



^en^bitiM) = V \Vlc?~9l {tiM) ■ 



L 



(65) 



/3/2 

7 ' 



G" (h,t) Ef (t, t2)+G> (h , t) (t, h 



G^ih, -it) i~iT,t2)+ G^ ih ,lT)^f (it, t2 ] 



ere ghiti^h) is the non-interacting Green's function 
of the uncontacted lead L, possibly subject to a uniform 
but time-dependent bias. The full Green's function in the 
central region will now obey the KBE with both the self 
energy from the interaction and the one from the leads 
• via the embedding. The embedding self energy is entirely 
non-local in time and will be treated on the same footing 
the correlation part of the interaction self energy. It 

rru 11- • • i 1 -4-1, f ii, i- 4- is worth noting that the embedding self energy involve 

Ihe collision integrals with one of the time arguments ° 6 sj- 

complex specialize to self-consistency and can thus be calculated once the 

external bias (if present) is known. 



/< {t, -ir) = [ diSf (i, t) G< {t, -ir) 
Jo 

h/dT[l]<(i,-zT) G^^(r-r) + E>(i,ir) G^^(-(T + r))] , 
Jo 



E. The Dyson equation for T and W 



I> i-iT,t) 



/3/2 



Jo 



dT [G^' (r - r) I]> (-ir, t) + G*^ (r + r) I]< (zr, t)] . 

(61) 



(60) Propagating in time the KBE within a specific MBA- 
based partial summation, requires solving a Dyson equa- 
tion for auxiliary quantities which enter the expression 
for the self energy. For the TMA and GWA such quan- 
tities are the T and W , respectively. The components of 
the corresponding Dyson equations (again specialized to 
the case of T) Eqs. ([STI [5^ are, for both times on the 
real axis. 



It is worth noting that all Eqs. ([551 - 1501) contain terms 
which involve integration along the Matsubara (vertical) 
axis and which represent the memory of initial state cor- 
relations during the time evolution. 

For the collision integrals, one can derive a similar sym- 
metry property: 



dt 



$^ (t 1 , t) UT^ {t, t2 ) + (i 1 , KT^ {t, t2 ) 



li/2 



/f (l2) = -/|(2l)^ 



(62) 



An important consequence of this relation is that the 
densities^ are manifestly real. From the KBE one can 
derive that the condition for real densities is given by 

[/< (t, t) - I< (t, t)] = [/< {t, t) - I< {t, t)] ^ , (63) 



- - / df[$< (ti ,-zf ) UT> {-if, ta) + ih , if) UT< (^f, ts)] , 

(66) 



(<1,<2) (^1,^2) 



-/ dt 
Jo 

r/3/2 



T^(ti, f)Zi$^(f, t2) + T${h,t) U^^{t, t2) 



which is manifestly satisfied by Eq. ((62)) . Furthermore, 
on the time diagonal, we obtain another relation, which 
is very useful in from the computationally point of view: 



1 ' 

- df[T<{ti-if)U^>{-if,t2)+T>{ti,if)U^<{ifM)\ , 

(67) 



^1/2 



^1/2- 



(64) when one of the time arguments is imaginary we have 



T< {t,-iT) = $< {t,-iT) 



dt<^'^{t,t)UT<{t,-iT) 



From the properties of integral equations, it follows that 
all the symmetry and structure properties of the Green's 
function are those of Gq. 

The discussion of the KBE above is valid for extended _ / df[^^ {t,—if)UT^^ {f — t) + {t,if)UT'^^ {— (f + t))] , 
and finite systems alike. When the interaction is con- Jo 

fined to a central region which is contacted to possibly (68) 
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T> i-iT,t) = $> {-iT,t) dtT>{-iT,f)U^^{t,t) 

Jo 

/•/3/2 

- / df [T^' (r - f) U<i>> (-zf , i) +r*^(T + f) W$< (if, t)] . 
Jo 

(69) 

F. Time propagation algorithm 



Up to this point we have not made any approximations 
but merely formal rewritings of the KBE. To calculate 
the S {Tp + A,Tp) we divide the interval [Tp,Tp + A] into 
N intervals in which the single-particle Hamiltonian, (), 
is taken constant and evaluated at the midpoint. For a 
constant t) we obtain 



As mentioned in section fVl Dl we treat the time-local 
part of the self energy {T,hf) on the same footing as 
the non-interacting terms (h) in the time propagation. 
The evolution from these single-particle terms can be ex- 
pressed in terms of a single-particle evolution operator 
S which is a time-dependent matrix in single-particle la- 
bels. This leads to the following unitary gauge transfor- 
mation 

{h,t2) = S (ti , 0) {h ,t2)S^ {t2 , 0) , (70) 

where S satisfies the following differential equation: 

idt,Sih,0) = l)iti)S{h,0), (71) 

with the initial condition 

5 (0, 0) S*^ (0, 0) = 1 (72) 

and the group property 

S ih,t) S it,t2) ^ S {h,h) . (73) 

Specializing to the case of G^, we get: 

^dt,G> {h,t2) - i) ih) S (ti, 0) g> {hM) ih, 0) 

+S{h,0)idt,g> {ti,t2)S^ (^2,0) 

^^{ti)G> {tiM) + Ii ihM)- (74) 

where the second equality comes from the Kadanoff- 
Baym equation for G^ . This results in 

idt,g> [hM) = SHti.Q) I> ihM) S {t2,Q) . (75) 

Therefore, by integrating from present time Tp to Tp-|- A, 

i[g> {Tp + A,t2)-g> {Tp,t2)] 

Tp+A 

{t,0)l> {lt2) Sit2,0)dt 

p 

/•A _ 

= {Tp, 0) S^t + Tp, Tp) I> {t + Tp, t2) S (i2, 0) cTt. 
Jo 

(76) 

Thus 

G> {Tp + A, t2) ^S{Tp + A, Tp) G> {Tp, ta) 

-iS {Tp + A, Tp) S^t + Tp,Tp) I> {t + Tp, ^2). (77) 
Jo 



(j + l)A J, , ^ 



N 



which is evaluated by diagonalization. The resulting ex- 
pression becomes 



S{Tp + A,Tp) ^Y[exp[-ii)[Tp- 



Given that A'^ is taken large enough, the only er- 
ror of the above expression comes from the extrapo- 
lation/interpolation of f), which is small as the den- 
sity (which enters f) via the HF term) is a continu- 
ous and smooth function of time. To solve Eq. ([77)) 
we also need to approximate the integral. This can 
be done in two different ways depending on which of 
the two quantities I^(i + Tp,t2) or I^(i + Tp,t2) = 
S'' {t + Tp,Tp) {t + Tp,t2) is the most slowly vary- 
ing function. We have tried both and seen that the 
!(" {t + Tp, ^2) is the smoothest. The integral is done nu- 
merically, typically with a 2- or 4-point formula. Similar 
expressions are used for the other KBE. Special atten- 
tion is needed only for the time diagonal. In this case we 
combine the two first KBE, Eqs. (|54l [SS]) and using the 
property in Eq. (|64|) . 



t[dt,+dt,]G< {h,t2) = [fi,G< {h,t2)] 

{h,t2) - If {h,t2) , (80) 

we then change to the variables t = (ii + t2) /2 and t' — 
ti — t2- This gives 

idtG<{t + t' I2,t~t' /2) = [f),G< {t + t' /2,t-t' /2)] 

If {t + t'/2, t - t'/2) - If {t + t' /2, t - t'/2){8l) 

By performing the same gauge transformation as above 
and setting i' = ti — t2 = we obtain^ 

idtg< {t) ^ i< {t) - i< (t) , (82) 

where /< (t) = S'f {t, 0) /< {t, t) S {t, 0). Integrating from 
time Tp to Tp -t- A we obtain 

g< {Tp + A)^g< {Tp) - z T"''^ (/< {t) - If (t)) di, 

(83) 

which leads to 

G< {Tp + A) = S{Tp + A, Tp) G< (Tp) {Tp + A, Tp) 

-iS{Tp + A, Tp) J (^if{t + Tp) - if{i+ Tp)) di 

xSHTp + A,Tp). (84) 



10 



The integral is then evaluated in the same way as dis- 
cussed above for the case ti ^ti. 



VII. TDDFT EXCHANGE-CORRELATION 
POTENTIAL FROM MBPT 

Our time-dependent densities from the different Many- 
Body Approximations also provide insight for the 
TDDFT exchange-correlation potentials for strongly cor- 
related systems. Given a specific approximation, from 
the resulting time-dependent density we obtain the cor- 
responding effective potential v^ff = Yj^i^w+Vxc^ where 
Efl- is the Hartree potential and Vxc the exchange corre- 
lation potential. In practice, this is done via a numerical 
reverse engineering procedurci^ This algorithm imposes 



a) n=l/2 



b) '« = //() '''''' 



that |n (i?, t) 



KS 



{R,t) 



n 



at each time-step, where 
lipl^^ {R,t) |2 is the Kohn-Sham den- 
sity and tp^^ are the Kohn-Sham orbitals. The n^^ is 
found by solving iip^^ = {i + v^ff) where the ki- 
netic energy is given hy t ~ —V'}2 



(Hi?': 



1 ' 1 ' 1 ' 1 ' 1 

a) n=l/2 n ,i 








___yvAA\GWA_; 


- ywv^ 






^__/AAa^_tma_^ 




ywW EX 





FIG. 4: (Color online) Ground-state spectral functions for 
Af=6 for weak interaction strength, (7=1, for different fill- 
ings. The curves correspond to exact (black), TMA (red), 
BA (green), GWA (blue) and HFA (orange). The curves are 
shifted for clearer comparison and we have broadened the dis- 
crete spectra with a Lorentzian with a FWHM = 0.4. 



VIII. RESULTS: GROUND STATE 

To start the time propagation of the KBE, one 
needs the initial, equilibrium, one-particle propagators. 
These are found by by solving the Dyson equation self- 
consistently. Here, we show only results for isolated clus- 
ters, and defer to a future publication the case of clus- 
ters contacted to leads. To characterize the ground-state 
properties, we present in Figs. |3]and[5]the spectral func- 
tions at the first, R = 1, site in the cluster. This is the 



jM 



HFA 




-4 4 8„0246 



10 



FIG. 5: (Color online) Ground-state spectral functions for 
M=6 for strong interaction strength, ?7 = 4, at different fill- 
ings. The curves correspond to exact (black), TMA (red), 
BA (green), GWA (blue) and HFA (orange). The curves are 
shifted for clearer comparison and we have broadened the dis- 
crete spectra with a Lorentzian with a FWHM — 0.4. 



site at which the external perturbation is applied (see 
section HXl) . The results are for the different Many-Body 
Approximations, at two particle concentrations and in- 
teraction strengths. A first, general comment is that the 
performance of MBA:s generally worsens with increasing 
U. Some basic, generic features of the different Many- 
Body Approximations we employ to evolve in time our 
clusters are already revealed by the ground-state results. 
Such features, specifically discussed for the cases in Figs, 
m and [51 are common to other clusters with different fill- 
ings, sizes and interactions strengths (such additional re- 
sults are not shown). Fig. prefers to the situation of 
weak interaction, where the performance of the MBA:s 
is satisfactory (in particular, all the MBA:s reproduce the 
symmetry breaking at low filling of the exact solution) . 

At half filling, a correlation gap opens on increasing 
U, which becomes clearly visible in the strong correla- 
tion regime. Fig. [S] This behavior in clusters is consis- 
tent with the exact solution of the ID Hubbard model in 
the thermodynamical limit4^ This feature is reproduced 
by the different approximate schemes. However, the size 
of the approximate gap depends on the Many-Body Ap- 
proximation: the best estimate is given by the BA value; 
yet, the latter largely underestimates the exact value. 
Missing the contribution of correlation effects, the non- 
magnetic HFA solutions reproduce the non-interacting 
spectral density. As a second general feature, the MBA:s 
introduce spurious satellites away from the band region. 
This problem is most pronounced in the BA and GWA 
curves. In the low density regime, for large U, a satellite 
structure appears in the exact solution (about 6.5 in the 
bottom of panel 46). In an extended system, this high 
energy spectral feature represents a two-electron anti- 
bound state (i.e. outside the band continuum). The 
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satellite is well reproduced by the TMA (although its 
distance from the band region is overestimated) while is 
smeared out in the BA and GWA and obviously absent 
in the HFA (as it includes no correlation effects). For the 
strong interaction case, the agreement of the different ap- 
proximations with the exact curve in the band region is 
only moderate. 

When comparing the SGWA to the GWA we see a 
slight improvement, which is expected, as the SGWA 
includes fewer faulty diagrams. This improved version 
of the GWA is, however, still worse than the BA or 
the TMA. In this section, we show no results for the 
SGWA, since it introduces only a a marginal improve- 
ment in the ground-state spectral density: however, we 
will present below SGWA time-dependent densities. It is 
worth noticing that, using a MBA expansion in terms of 
non-magnetic propagators, the SGWA will have a mag- 
netic instability on increasing U. This can be seen most 
easily in a Hubbard dimer with two electrons with oppo- 
site spins. In the dimer, where the poles of W^[G'o] are 
e = ±V4F^ ± 2VU, this unphysical symmetry break- 
ing occurs for U > 2V. Conversely, the exact ground 
state for a dimer is always a spin S = (singlet) state, 
since, for any positive U, the dimer ground-state energy 

gs ^ dimer 2 U ^ dimer' 

It is worth noting that all the spectral functions shown 
in Figs, m and [5] are not as good as those obtained with- 
out self-consistency, i.e. when stopping after the first 
iteration of the Dyson equation, see Fig [6] This is an ex- 
ample of the known fact that self-consistent conserving 
approximations often have worse spectral properties than 
non self-consistent ones J^iii To guarantee the fulfillment 
of the conservation laws and to get unambiguous total 
energy results it is capital to achieve self-consistency. 
If the quantity of interest is the spectral functions one 
should instead make use of different partial summation 
criteria)^ and in some cases include vertex corrections 
to remove artifacts introduced by self-consistencyi^ 

As a final remark to this section we note that all MBA:s 
based on partial summations involve infinitely many pos- 
sible excitations. These excitations, represented by dia- 
grams in the self energy, result in infinitely many, but 
discrete, number of poles in the ground-state spectral 
functions4^ The exact solution, in contrast, lives in a 
finite phase space which implies a finite number of poles. 



IX. RESULTS: TIME DEPENDENCE 
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FIG. 6: (Color online) First iteration and self-consistent TMA 
spectral function versus the exact one. M=6 and in a): n — 
1/6 and in 6): n = 1/2. The curves correspond to exact 
(black), first iteration TMA (green) and self-consistent TMA 
(thick blue). The curves are shifted for clearer comparison 
and we have broadened the discrete spectra with a Lorentzian 
with a FWHM = 0.4. 

existence of multiple steady states in isolated and con- 
tacted clusters. 

Isolated clusters: MBA:s vs exact results. We start the 
time evolution at t = with the ground-state Green's 
function. For positive times t > we apply a spin- 
independent external field to the system. We have stud- 
ied different types of external fields, but in this paper we 
present results only for the form wr (t) — woSii,i&{t), 
that is we consider a step perturbation and let it to act 
only on the leftmost, R = I, site. The time is given 
in units of the inverse hopping parameter (l/V) and 
all curves represent the dynamics on site R = 1. The 
cases displayed are the same as those considered for the 
ground-state results. 
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In this section we examine the performance of the dif- 
ferent MBA:s. To accomplish this, we use as benchmark 
exact many-body solutions. In general, the latter are 
available only for finite systems, and numerical in na- 
ture. As a consequence, in this section we deal exclusively 
with isolated clusters, and focus on general aspects of the 
time-dependent densities and MBA:s. However, we defer 
to the next section two important outcomes of the KBE 
time evolution: the correlation-induced damping and the 



FIG. 7: (Color online) Time-dependent densities on site 1 for 
M = 6. The curves correspond to exact (thick solid black), 
TMA (dashed red), BA (dotted green), GWA (thin solid blue) 
and HFA (brown dashed dot). 

We show the resulting time-dependent densities in 
Fig. [T] The curves correspond to some initial states 
shown in Figs. |3]and[5] In the panels a) and b), n = 1/2 
while in panels c) and d), n = 1/6. In the simplest case 
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{U = l,wo = 1), presented in panel a), all MBA:s give 
a good description of the density. On increasing the 
strength of the interaction and the external field, panel 

b) , we clearly see that the HFA description is rather 
crude, whilst the curves from the other MBA:s are very 
similar to each other and closer to the exact density. 
In the case of strong interaction but weak field, panel 

c) , we see that none of the MBA:s give an adequate 
description. We interpret these results as a consequence 
to the fact that the ground-state spectral function is 
not well described in the band region: the latter is 
responsible for the response to weak fields. On the 
contrary, when the field is strong, panel d), the TMA 
performs much better than the other MBA:s. This is 
due to the fact that the non-linear response involves 
states at higher excitation energy and the TMA is the 
only approximation which, to some extent, reproduces 
the satellite structure. 




FIG. 8: (Color online) Time-dependent Veff on site 1 for 
Af = 6, (7 = 4 and wo = 5. The curves correspond to exact 
(thick solid black), TMA (dashed red), BA (dotted green), 
GWA (thin solid blue) and HFA (brown dashed dot). 

Isolated clusters: MBA:s, exact results and TDDFT. 
It is interesting to examine some of the results just pre- 
sented from a TDDFT perspective. A clear advantage 
of TDDFT is that, for the time evolution, it deals with 
quantities with a single time argument (one is propa- 
gating the Kohn-Sham orbitals). Nevertheless, a key re- 
quirement in TDDFT is that v^c (an thus We//) should 
depend in a non-local (in space and time) fashion on the 
particle density. In this way, all the complexities of the 
many-body dynamics arc subsumed into a highly non- 
trivial dependence of Vxc on the density. Not surpris- 
ingly, making progress in the construction of improved 
XC functional is a rather challenging task, especially for 
strongly correlated systems. It is also true that, in some 
cases, simple adiabatic local approximations can provide 
satisfactory results; this is also the case for a TDDFT de- 
scription of the Hubbard model in non equilibrium.— 
However, memory and non-local effects should in gen- 
eral be taken into account. The so-called variational ap- 
proach to TDDFT^i^ has the advantage of a systematic 



inclusion of many-body contributions in the XC poten- 
tial. Also, in this way non-locality in space and memory 
effects are properly included, once Vxc is retrieved from 
the many-body self energy via the time-dependent Sham- 
Schliiter equation4^ Deferring a "bottom- up" construc- 
tion of Vxc via MBA:s to future work, we here still wish 
to explore the connection between TDDFT and MBA:s 
on the Keldysh contour, looking at exchange-correlation 
potentials obtained via time-dependent reverse engineer- 
ing using the time-dependent densities from the KBE. 
The results of this procedure are presented in in Fig. O 
M = 6, N = 2,U = 4, Wo = 5. This corresponds to the 
time-dependent densities presented in panel d) of Fig. [T] 
Consistently with the density results, we see the v^ff in 
the TMA is superior to those from the other MBA:s, and 
quite close to the exact one. We also note the large dis- 
crepancy of the HFA and that, for both the BA and the 
GWA, Veff exhibits a damped behavior (the same is ob- 
served in the densities in panel d) of Fig. [71 see section 
|X] below). In spite of not being perfect, the agreement of 
Veff from the TMA with the exact one is quite encour- 
aging, suggesting that there is ample scope for pursuing 
the construction of improved v^c from (suitably chosen) 
MBA:s. 




FIG. 9: (Color online) Time-dependent densities for M = 
4, U — 1.5, lOo = 5 and n — 1/4. The curves correspond 
to exact (black), GWA (dotted red), SGWA (thick blue) and 
TMA (dashed green). 

Spin-dependent GWA. Before concluding this section, 
we wish to discuss briefly the effect on making a spin- 
dependent treatment in the GWA. In Fig. [5J we make 
a comparison between GWA. SGWA and, for reference, 
TMA. Due to the symmetry breaking of the SGWA (see 
section rVIIip . and that we only consider the spin unpo- 
larized case (otherwise, one should start with polarized 
propagators) we confine ourselves to the weak interac- 
tion regime. As evident from Fig. [51 for the chosen 
parameters the SGWA is slightly better than its spin- 
independent counterpart but inferior to the TMA. 

As a final remark to this section, we note that in the 
general case, GWA was designed to give a good screen- 
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ing of the Coloumb interaction^ which in our system has 
aheady been taken into account indirectly by the model 
itself. The TMA, on the other hand, is known to give a 
good performance if the interaction is short ranged, es- 
pecially in the low density regimei^ The general good 
description of the the TMA (both in and out of equi- 
librium) for the short-ranged Hubbard Hamiltonian is in 
accordance to previous studies of ground state properties 
of clustersi^i^ One can in other words say that the per- 
formance of the different approximations in the ground 
state has non negligible relevance to the time-dependent 
behavior of the system, and this is especially the case 
for two most distinctive spectral features, that is: the 
band-gap and satellite structure. 



X. RESULTS: DAMPING AND MULTIPLE 
STEADY STATES 




FIG. 10: (Color online) Densities for M = 2, n = 1/2, i7 = 
1, Mo = 5, exact (black), GWA (dashed red) and SGWA 
(blue) . In a) : damping of GWA density versus exact solution. 
In 6): time-dependent densities for the GWA, initialized with 
the self-consistent GWA ground state and the non-interacting 
Go. 

In this section we will investigate the correlation- 
induced damping and the multiple solutions of the sta- 
tionary KBE. We find that these features are general for 
all MBA:s which include correlation effects. We exem- 
plify this by presenting results for different MBA:s in the 
various parts of this section. 

When we let our isolated systcm(s) evolve under the 
action of a strong external field we reach an artificial 
steady statej^i^ see Fig. rTOf a). The damping mecha- 
nism is not a mere consequence of the infinite number of 
poles in the initial state but is rather intrinsically linked 
to the time propagation scheme. To show this fact we 
present in Fig. \TU[b) the time evolved density initiated 
with the non- interacting propagator which has a 

finite number of poles. ^From the curves in Fig. llOf b') it 
is evident that the non-interacting initial state leads to a 



very similar damped density profile. Note that the sim- 
ilarity of the curves in Fig. [TUr 6). indicate a robustness 
of the KBE time evolution against the initial conditions. 




FIG. 11: (Color online) Conservation laws and time reversal 
in the GWA for M = 2, n = 1/2, [/ = 1, wo = 5. 

The damping is not a numerical artifact: in Fig. [TT] 
we see that particle and energy conservation are strictly 
obeyed within our numerical accuracy, see top and mid- 
dle panel. Moreover, the evolution satisfies time-reversal 
symmetry. That is, when we reverse the direction of time 
in the propagation, the system goes back to the initial 
state and remains there, see bottom panel. The damp- 
ing rate increases with the strength of the external field 
and is absent in the regime of linear response. The dy- 
namics in this limit is described by the Bethe-Salpeter 
equation, with a kernel 5T,/ 5G. The latter would have a 
discrete spectrum in our MBA:s, and so would the result- 
ing density response. A discrete response function will in 
turn lead to a non-damped dynamics. We wish to stress 
that in a formulation based on conserving Many-Body 
Approximations, the key quantity is the generating func- 
tional; the corresponding G is then defined via the self- 
consistent KBE, and not via an underlying wave function 
scheme. Without the connection to the underlying wave 
function, there is no guarantee that, for example, systems 
with a finite phase space will only have a finite number 
of excited states or that they will not have a damped 
dynamics. 

To study the non-linear response regime, we find con- 
venient to introduce the instantaneous spectral function: 

A (Tp, w) = -Tr Im j e'"^ [G>- G<] (Tp+-, Tp--j dr, 

(85) 

where Tp ~ (ti + ^2) /2 and t ^ [ti — 12) , and its coun- 
terpart in time space. When reaching the steady state, 
the spectral function gets broadened in energy space. Fig. 
[T^g). and damped in time space. Fig. [T^6). Note that, 
for our dimer, the exact instantaneous spectral function 
would continue to oscillate in time space. The differ- 
ent MBA:s have different damping rates; among them, 
the TMA is in general the slowest. The damping acts 
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such perfect compensation. The infinite number of poles 
of the ground state and the correlation-induced damping 
will thus be artificial for a finite system. 



FIG. 12: (Color online) Time-dependent spectral functions in 
the GWA for M = 2, n = 1/2, U = I, wo = 2, a): in energy 
space and 6): in time space. The different curves correspond 
to Tp = 10 (black), Tp = 35 (red) and Tp = 100 (blue). 



strongest on the perturbed site, it generally decreases 
with system size and is most important at half filling. In 
our study we have not made an exhaustive characteriza- 
tion of how the large-size limit is gradually obtained; this 
is indeed an the issue we plan to address in future work. 
Still, we wish to present in the rest of this section some 
general considerations on the effect of system-size on the 
damping behavior. 

System size and damping. In an exact treatment, a fi- 
nite system has a corresponding finite phase space and 
will thus not be able to fully relax to a stationary steady 
state. In a large but finite phase space it will, due to deco- 
herence, give rise to a quasi steady state. The larger the 
system gets, the more and more complete the damping 
becomes. Thus, after a long time, in an exact time dy- 
namics the system would exhibit noise-like fiuctuations 
which never die but which decrease in amplitude with 
system size. We have seen that in our approximate KBE 
evolution the system, while being finite, attains a sta- 
tionary state. Thus, the artificial damping of the KBE 
dynamics is expected to increasingly resemble the physi- 
cal damping as the system grows. 

Self- consistency and damping. In Many-Body Perturba- 
tion Theory, the self energy accounts for possible exci- 
tations of the system which involve a certain number of 
particles/holes. Any MBA based upon partial summa- 
tions includes diagrams of all orders. These terms act 
as an effective bath which gives rise to infinitely many 
discrete poles in the ground state (in isolated systems) 
and correlation-induced damping in the time dynamics. 
In a finite system, there will be contributions to the self 
energy which annihilate more holes/particles than those 
which can be accommodated in the system. In an exact 
theory, these unphysical terms would be exactly canceled, 
order by order, by other unphysical pieces. In approxima- 
tions such as the GWA or the TMA, there is in general no 




FIG. 13: (Color online) Time-dependent densities for Born 
approximations with different levels of seff-consistency for 
M = 2, n = 1/2, U = I, Wo = 5; BAo (black), BAhfa (red) 
and BA (blue). 

In Fig. 1131 we present results from three versions of 
the BA to illustrate the effect of an increasing level of 
self-consistency. These approximations, which arc all 
particle-conserving)^ involve different polarization prop- 
agators. In the first case we evaluate the polarization 
propagator with ground-state propagators (BAq) (top 
panel). In this case the density does not damp. In the 
second case we evaluate the polarization with propaga- 
tors in the time-dependent HFA approximation (BAjjfa) 
(middle panel). In this case we observe partial damp- 
ing. If we finally use the self-consistent G (BA) (bottom 
panel), we get complete damping. We see in other words 
that if all G:s that build up the the self energy are the 
self-consistent ones we reach a steady state.— 
Multiple .steady states. Another striking feature related 
to the correlation-induced damping is that the steady 
state is not unique for a given final external field: it de- 
pends on how the perturbation is switched on. In our 
simulations, in the case of adiabatic turn-on, we reach 
the ground state of a system with an on-site energy cor- 
responding to the final external perturbation. This is 
consistent with the adiabatic theorem. If, however, the 
perturbation is switched on suddenly, we reach a non- 
physical steady state with the same energy as at i = 0+ . 
This non uniqueness is indicative of an important as- 
pect: given a final external potential, there are multiple, 
in principle infinitely many, solutions of the stationary 
KBE in finite clusters. Furthermore, as we show next, 
this general statement, in some cases, remains valid in 
presence of macroscopic leads, where we have both a self 
energy Ttemb from the leads and a self energy "Emba from 
the interactions. Both contributions are non-local in time 
and may lead to damping. Thus, in general we have two 
damping mechanisms: One due to the coupling to the 



15 



0.5 




10 20 30 40 

time 



FIG. 14: (Color online) Time-dependent densities for a cluster 
isolated and coupled to leads in the TMA with sudden and 
slow switch-on. M = 2, U = 1, wo — 8. Lead parameters: Vl 
= 1, Vlc = 0.5, filling = 0.5 and bias = 0. a) Sudden switch- 
on: isolated cluster (thin black), coupled to leads (thick red). 
b) Slow switch-on: isolated cluster (thin black), coupled to 
leads (thick red); ground-state densities of final Hamiltonian: 
isolated cluster (dashed black), coupled to leads (dashed thick 
red). Inset: expanded view of the long time limit behavior. 



continuous lead band and one induced by correlations. 

In Fig. [14] we present the time-dependent densities 
within the TMA for a dimer, both isolated and coupled 
to unbiased leads f^i^ subject to an external field with 
sudden and slow {w (t) = wot/tmax) switch-on, (panel a) 
and panel b) respectively). 

The strength of the perturbation is such that the 
ground state corresponding to the final external poten- 
tial contains states outside the band continuum (bound 
states); in this case we find that the KBE admit mul- 
tiple steady-state solutions<^ These steady states vary 
continuosly with the way the final potential is reached. 
In other words we have infinitely many steady states 
where similiar switch-ons give similar steady states. In 
the case of the HFA we see that the density continues 
to oscillate if we have more than one pole in the ground 
state of the final Hamiltonian. This results is consistent 
with recent work on the role of bound states in quantum 
transport 1^^— However, we wish to remark that, in the 
HFA, the frequency, amplitude and average value of the 
oscillating density will, however, depend on the way the 
external field is switched on<^ 

Similarly to the isolated case, we find that for a slow 
switch-on we tend to reach the ground state correspond- 
ing to the final Hamiltonian (given by the dashed curves), 
while for the sudden switch-on we reach another steady 
state. If the external field is such that the corresponding 
ground state of the final Hamiltonian includes only exci- 
tations within the band then the steady state reached is 
the final ground state, independently on how the pertur- 
bation is switched on. 

The damping induced by correlations and the one due 




time 



FIG. 15: (Color online) Time-dependent densities for a cluster 
with different couplings to unbiased leads in the BA. M = 2, 
U=1,wo^5,Vl= 5. 



to the coupling to the macroscopic contacts have in gen- 
eral different characteristic time scales. This is clearly 
seen in Fig. [15] where we have gradually increased the 
coupling to the leads: here, the external field and the 
lead band width are such that the corresponding ground 
state of the final Hamiltonian contains no poles. In the 
case of an isolated cluster {Vlc = 0), the damping is 
only due to correlations. Once the coupling to the leads 
is non-zero, there will be, in addition to the correlation- 
induced one, a damping due to the contacts; the lat- 
ter will eventually bring the system to the corresponding 
ground state. In Fig. [15] , the ground-state values of 
the densities are represented by arrows (on the scale of 
the figure, the red and green arrows are not distinguish- 
able). Note that each curve in Fig. [TSl corresponds to 
a different system, i.e. with a different device-lead cou- 
pling strength. Each of these systems has a different final 
Hamiltonian, i.e. a different corresponding ground state. 
When the couphng strength is weak {Vlc = 0.125,0.5), 
the time scale of the damping due to the leads is much 
longer than correlation-induced one. For intermediate 
couphng strengths {Vlc = 1-5), the characteristic times 
will be of the same order of magnitude and the interplay 
of the two mechanisms becomes intrisically difficult to 
discern. When the coupling strength is strong {Vlc = 3), 
the damping is completely dominated by the leads. 

As a final remark, we know that the correlation- 
induced damping is artificial in isolated clusters, and we 
saw in Fig. [15] that when the coupling to the leads is 
weak, the initial damping is also dominated by correla- 
tions. These two facts together seem to cast some doubt 
on the capability of the KBE-|-MBA:s scheme to describe 
(within the simple MBA:s discusses here) transients in 
the weak coupling case (a regime which is of high experi- 
mental interest, and in fact among the most investigated 
in the literature). 

To summarize this section, we see that correlation- 
induced damping and multiple steady states are present 



16 



both in isolated and contacted clusters. In the case of 
isolated clusters, this damping is artificial as the exact 
solution does not reach a steady state. The question 
whether or not the correlation-induced damping and the 
existence of multiple steady states in a cluster coupled to 
leads is an artifact of Many-Body Perturbation Theory is 
beyond the scope of this paper and is left to future work. 



XI. CONCLUSIONS AND OUTLOOK 

A main objective of this paper has been to describe 
in detail a method, within the framework of the time- 
dependent Kadanoff-Baym equations (KBE), to study 
finite systems in equilibrium as well as out of equilib- 
rium. The main emphasis of the paper has been on finite 
clusters, for which a meromorphic representation of the 
equilibrium many-body quantities is possible, but, in few 
instances, we have also considered some results for con- 
tacted clusters. 

As a concrete application, we examined the time evolu- 
tion of clusters with strong, short ranged electron inter- 
actions within several Many-Body Approximations, and 
compared the results to exact ones. A first main outcome 
of this comparisons is that, for short ranged, Hubbard- 
type interactions, the T-matrix approximation performs 
very well at low densities, and is in general superior to the 
GW and second Born approximations, both in describ- 
ing the time-dependent density and the corresponding 
exchange-correlation potential. 

A second important outcome of our work is the exis- 
tence of two remarkable features of the time-dependent 
KBE. The first is that the KBE present a correlation- 
induced damping in the non-linear regime. The second is 
that the steady state reached is in general not unique, i.e. 
the stationary KBE support multiple stationary states. 
Since a finite cluster subject to a non-adiabatic pertur- 
bation will oscillate indefinitely, it is clear that, for finite 
systems, such damping and multiple steady-state behav- 
ior are artificial. We argue that this shortcomings will 
always be present when applying infinite order pertur- 
bation theory based on partial summations to finite sys- 
tems. In this paper, we were not able to provide a con- 
clusive answer to whether such two aspects of the KBE 
bear any physical meaning or if they are a spurious effect 
of MBPT. As we discuss next, this is part of our plans 
for follow-up work. 

Future research activity may be envisaged in various 
directions. On the methodological side, it would be of 
interest to further investigate, for contacted systems, how 
sound the correlation-induced damping and the multiple 
steady states arei^ 

A second methodological issue we are currently ad- 
dressing is performance of the different MBA:s in quan- 
tum transport geometries. We will will do by compar- 
ing the approximate results with those of time-dependent 
density-matrix renormalization group calculationsi^ 

Another possible line of study would be searching for 



algorithms to reduce/optimize the computational costs 
of time-dependent KBE numerical calculations. This is 
necessary to deal with more realistic systems. 

As a fourth direction, it would be of great interest to 
study the effect of including bosonic degrees of freedom as 
vibration phenomena often play an crucial role in quan- 
tum transport. 

Finally, we also intend to use our KBE-based compu- 
tational treatments in some specific applications: Pos- 
sible examples are real-time qubit manipulation, bista- 
bility induced by electron-phonon interaction, cold-atom 
dynamics. More specific details are deferred to future 
publications. 
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APPENDIX A: Energy functionals 

The energy functionals considered here are the 
Galitskii-Migdal (GM) and the Luttinger-Ward (LW). 
An important difference between them is that the GM 
needs only the knowledge of the single-particle Hamil- 
tonian and of the spectral function of G, while the LW 
incorporates the $ functional as well as S [G] , which de- 
pends on which approximation one is using. Here we 
give the expression for the LW functional for the GW- 
approximation. In practice, we evaluated the LW only in 
the ground state, i.e. in frequency space, while the GM 
was also used during time evolution. 

a. Galitskii-Migdal 

In frequency space, the GM functional reads^ 

E = -'-Ti-{{e + h)G}, (A-1) 

where Tr stands for traced and h is the single-particle 
Hamiltonian not including T,hf- Thus 

riF ^F 

where np is the number of poles below the Fermi energy, 
and after having made use of the meromorphic represen- 
tation of . In time space, the form of E becomes 

E = -'-Tr [{idt + h) G< {t, t+)] . (A-3) 
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Making use of the the equation of for G, we obtain 
E = -'-Tr[hG< {t,t+)] 

-^Tr [{h + ^hf) G< (t, t+) + /< (t, . (A-4) 



b. hutting er- Ward 

The LW energy funetional is^ 

= $ - Tr {SG + In (S - Gq ^) } , (A-5) 

where the self energy depends functionally on the input 
G i.e. S = E [G] and $ is the generating (of the MBA:s) 
functional. Specializing to the case of the GW approxi- 
mation, the $ functional becomes 

^GW = $HF + ^Tr {In [1 - UP] + UP} , (A-6) 

where ^hf = («/2) Tr {S/f^iG] G}. When computing 
the logarithmic term In (S — G(^^) in Eq. (jA-5[) . it is 
useful to make the following separation 

In [S - Go'] = In [-G^'p] + In [1 - Ghf^c] , (A-7) 



where Gj^p = e — (i and f) incorporates the the Hartree 
and Fock terms and where Sc is the correlation part 
of the GW self energy. The first term in Eq. (|A-7P 
is the sum of the eigenvalues of the occupied Hartree- 
Fock single-particle states (calculated from the correlated 
Green's function) which obviously does not correspond to 
the Hartrce-Fock energy. The second contribution in Eq. 
(|A-7P can be evaluated with a well known identitj*^: 



lull -Ghf^c]^- 



Ghf'^c 



1 — XGhf^c 



G^rdX. 



-dX 



(A- 



In Eq. ()A-8p . G satisfies the Dyson equation 



G = Ghf + XGhf'^cG. 



(A-9) 



The other logarithmic term in $gw, i- c. In [1 — UP], 
is treated in a similar way. In the actual calculation, ex- 
pressions of the form Tr [AB] are evaluated analytically, 
whilst integrals of the kind ABdX are performed nu- 
merically. 
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